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P H E b' A C bC 


The material included in this internal note was developed 
by William M. Lear of TRW under contract with the Mathematical 
Physics branch (MPB) at the Johnson Space Center (JSC). The 
material was presented to MPb in the form of a TRW publication, 
number 30759-6003-TU-00 , dated June 1, 1 977. The material was 
prepared for JSC publication by Paul H. Mitchell of the MPb, 
and the text is taken from the referenced TRW document in its 
entirety, with the following changes. 

(a) The coefficient of n in equation 13 (page 4) was 

changed from (C 22 + 2 C 12 ) (^22 consistent with 

the expansion of equation 12 and the following definitions. 

(b) First person pronouns were changed to third person 
pronouns on pages 11 and 18. 

(c) The concluding remarks were labeled as "8. CONCLUSIONS" 
on page 19. 
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LEAST SQUARES POLYNOMIAL FITS AND THEIR ACCURACY 
By William M. Lear 
TRW Systems Group 
1 . INTRODUCTION 


Suppose we are given a table of data as 
shown at the right. And, suppose we desire 
to make a least-squares polynomial fit to 
this data. We will assume that t^-t^ ^ = aT 
Is a constant. That Is, we have evenly 
spaced tabular values of y*. Note that 


n 

t 


1 


1 

A 1 

2 

‘2 

A 

3 

^3 

A 



1 

1 

. 


• 

• 


4 

N 

^N 

A 


We will assume that n Is the Independent variable In what follows. To make t 
the Independent variable we use Eq. 1 to substitute for n. 

For the sake of explanation, let us assume that y* Is given by the second 
order polynomial 


yj - do + d,n + (2) 

where q Is a zero-mean, timewise uncorrelated random variable (noice) whose 
n 2 

variance is o^. Independent of n. In order to predict future measurements, 
we need to know a^, a-|, and a 2 Thus our state vector Is 


X 



(3) 


Using ~ to Indicate estimated value, we have for our estimated measurement 


(4) 


1 


The measurement partial derivative matrix, P = will be 


P = 1 


( 5 ) 


1 N N J 
Note that we have a linear measurements problem, 

£ = Px 


( 6 ) 


That is, the estimated measurement vector is a linear combination of the 
state vector elements. 

The least-squares solution for £ is given by 


(P^P)' 


IpT 


(7) 


X = £ + W(y^* - (8) 

C = Oq(P^P)‘^ (9) 

where Eq. 7 is the equation for the residual weighting matrix, and Eq. 9 
is the equation for*the state error covariance matrix, C = E [(£-£) (£-£)^] . 

Eq. 8 for x is iterated until converged. But from Eq. 6 and Eq. 7, Eq. 8 

becomes 


2 



or 



V 


r. 


X = X + - Px) 

= X + (P^P)"VV - (P^P)’V^Px 


X 


x = (P^P)’^pV 


( 10 ) 


which requires no iteration to find and is the least-squares equation 
for solving for the state vector when £ = Px. In the following sections we 
will give the analytic expressions for (P^P)"^ for various orders of polynom- 
ials. 

Now y* is given by 



+ q 


n 


= a„ + a,n 
0 1 


+ a2n 


(11) 


where y is the error free value of the ineasurement. Suppose we wish to know 
n 

the accuracy of our estimate of y^. 

E[(yp-y^)^] » E[((V^o^ (a^-a^)n + (a2-a2)n^}^] (12) 

= E[(a^-a^)2 ^ (a^-a^)V + (52-a2)^n^ 

+ 2(aQ-a^)(a^-a^)n + 2(aQ-a^)(a2-a2)n^ 

+ 2(a^-a^)(a2-a2)n^] 

But 

E[(V3o)^] = 

E[(a-| -a-| ) ] ~ C 22 
E[(®2"®2) ^ * ^33 


3 




^12 ^21 



^ ^13 ^ ^31 

E[(a^-ai )(a 2 - 32 )] = C23 = 0^2 

Thus we have 

' E[(y„-y„)^] = c,, + 2 C, 2 n t (C 22 + 2 C, 3 )n^ (13) 

t 2C2jn^ + Cjjn'* 

where., we remember that the matrix C was given by 

C = o^(p'^P)''' 

The C . . values will be functions of the total number of measurements, N, and 
2 ^ J 

a„. Eq. 13 will tell us how much we are able to reduce the measurement 

q ^ 2 

noise variance when we make a least-squares fit to obtain y^. We will be 
particularly interested in the values of E[(y^-y^)^] at n = 1, n = N, and the 
midpoint n = (N+1 )/2. 

We may also be interested in how accurately we can predict velocity and 
acceleration in our above example. Observe that from Eq. 1 


and 


From Eq. 12 


dy _ dy dn _ dy 1 
dt ~ dn dt ” dn aT 



. 2 
dy dn d y 1 


^ ' It (^i+^a^n) 


V 7 

dt^ aT^ 


(14) 

(15) 


I 

) 



Thus 


E[(yp"yj^)^] = E[{(a^-ai) + 2(32-32)0}^] 

at 


* ~^E[(a^-a^)^ + 4(82-32)^0^ 

at 


Or 


+ 4(a^-a^)(a2-a2)n] 


E[(yn-yn)^^ = ^ (C22 + 40230 + 40330^) 


aod likewise 


( 16 ) 


° ^ *^33 

lo the succeediog sectioos, the foil owl og equatloos have beeo used io 
the evaluatioo of P^P. 

1+2+3+...+N = ^ (N+1) 

1^+2W+...+N^ = ^N+1)(2N+1) 

l^+2^+3^+...+N^ = 5p(N+l)^ 

l^+2^+3^+...+N^ = ^ (N+l)(2N+i)(3N^+3N-l) 

2 

l^+2^+3^+...+N^ = ^ (N+1)^(2N^+2N-1) 

!®+2®+3®+...+N° = (N+1)(2N+1)(3N^+6N^-3N+1) 

l^+2W+...+N^ = ^ (N+l)^(3N^+6N^-N^-4N+2) 
r+2®+3”+...+N^ = (N+1)(2N+1)(5N®+15N®+5N^ 

-15N^-N^+9N-3) 


5 


V 


2. ZERO ORDER POLYNOMIAL, 


Our first example is for a zero order polynomial 


P = 


1 


(P^P)'^ = I 


(18) 


(19) 


E[(;„-y„)'] = 


( 20 ) 


3. FIRST ORDER POLYNOMIAL, >* * ♦ a,n 

n 0 I 


In this case 


P = 


1 1 
1 2 
1 3 




1 N 


( 21 ) 


(P^P)"^ 


N(N‘-1) 

2 

NCN'^-l) 


(N+1)(2N+1) -3(N+1) 

-3(N+1) 6 


2o* 

E[(y -yJ^] = [(N+l)(2N+l)-6(N+l)n+6n^] 

" " N(N-l) 


( 22 ) 


(23) 


2o 


= (2N-1) for n = 1 , n 


= Oq/N for n = (N+l)/2 


(24) 

(25) 


Note that n » (N+l)/2 yields the minitnum value for E[(y^-y^) ]. That is, 
our best value of y^ occurs at the midpoint of the fit interval. 


12o 


E[(y.-yJ ]= , , 

" ” aT‘'N(N*^-1) 


~ 2 ^2 
r^&i/ 


( 26 ) 


Notice that the total time, T, over which the fit Is made Is (N-l)aT, 


4. SECOND ORDER POLYNOMIAL, = a^+a^ 0+620^ 

1 1 1^1 

1 2 2^' 

P = i 1 3 3^ . 


1 N N 


( 27 ) 


(P^P)"^ 


N(N^-l)(N^-4) 


9 3o„ 9 

k ^ T ^ M r/fttin \ / 


(N+1)(N+2)(3N'^+3N+2) -6(N+l)(N+2)(2Nfl) 10(N+l)(N+2) 

-6(N+1)(N+2)(2N+1) 4(2N+1){8N+11) -60(N+1) 

10(N+l)(N+2) -60(N+1) 60 

(28) 

(29) 


E[{; -y J‘] = g ^ "2 [(N+1)(N+2)(3NS3N+2) 

" " N(N -l)(N^-4) 


-12(N+l)(N+2)(2N+l)n + 12(7N^+15N+7)n^ - 120(N+l)n^ + 60n^J 


3a' 


" N ' ( ' N+? ) ( N+2 ) -3N+2) for n = I , n = N 

2 


3o 


4N(N^-4) 


(3N'^-7) for n = (N+l)/2 


1 2a 


E[(y„-yn) ^ P 7 [(2N+l)(8N+ll)-60(N+l)n+60n^] 

" " ArN(N'^-l)(N^-4) 


12a 


—5 5 (16N'^-30N+11) for n = 1, n = N 

at N(N -1)(n'^-4) 


12a' 


—2 — for n = (N+l)/2 
at N(N -1) 


(30) 

(31) 

(32) 

(33) 

(34) 


8 



Note that Eq. 34 Is identical to its counterpart for the linear fit, Eq, 26. 

(35) 

It is of interest to note that the maximum accuracy of is not at 
the midpoint of the fit interval. This can best be seen for N >> 1 . Eq, 

29 becomes 



E[(y„-y„)‘] 


at N(N -1)(N'^-4) 


For J ^ = 0 E[(y^-y^)^] = ^ (a maximum) (37) 

For J ^ ^ E[(yn-yn)^] * ^ minimum) (38) 



where 


5. THIRD ORDER POLYNOMIAL, y = a +a,n+a,n^+aon^ 

n 0 1 d 3 

1 1 

1 2 2^ 2^ 

1 3 3^ 3^ 

P = 

I 

1 . . . 

2 3 

1 N N N'^ 


(39) 


(P'^P)’^ = 2 ^ 2 

N(N -1)(N'^-4)(N -9) 


11 

Ai2 

^13 

^4 

'21 

A22 

^23 

^24 

'31 

^32 

^33 

^34 

'41 

A42 

^43 

A44 


(40) 


= 8(N+l)(N+2)(N+3)(2N+l)(N^+N+3) 

Ai 2 = A2^ = -20(N+l)(N+2)(N+3)(6N^+6N+5) 
A^3 = A^^ = 120(N+1)(N+2)(N+3)(2N+1) 

^14 = -140(N+1)(N+2)(N+3) 

A22 = 200(6N^+27N^+42N^+30N+11) 

A23 = A32 = -300(N+l)(3N+2)(3N+5) 


X 


I 


t 


I 

* 



i 


10 


^24 " ^42 ' 280(6N^+15N+11) 
A 33 = 360(2K .1)(9N+13) 

A 34 = A ^3 = -4200(N+1) 

A 44 = 2800 


The accuracy of is given by 


E[(y„-yJ^] = 5 — [(N+l)(N+2)(Nt3)(2N+l)(N^N+3) 

-5(N+l)(N+2)(N+3)('5N^+6N+5)n+5(42N^+213N^+378N^+288N+91)n^ 

(41) 


-1 0( N+1 ) ( 71 N^+1 75N+96)n^+5( 246N^+525N+271 )n^ 


-1050(N+l)n^+350n®] 
8o' 


' (2N^-3N^*7N-3) for n • 1 . r - N (42) 

2 


3o 


4N(N‘^-4) 


A- (3N^-7) forn=(N+l)/2 

C A\ 


(43) 


Note that Eq. 43 is identical to its counterpart in the preceding section, 
Eq. 31. In this case however, it is suggested that we have a minimum value 
this time instead of a maximum value as before. 


11 


E[(y„-yn) ] 


40oJ/4T^ . , 2 

5 9^ 5-- [5(6N>27N'^+42NS30N+11) 

N(N -1 )(N^-4)(N -9) 


-30( N+1 ) ( 3N+2 ) ( 3N+5 )n+6{ 1 50N^+31 5N+1 55)n^ 
-1260(N+l)n^+630n^] 


For n = 1 and n = N 


iHl -y )^] = 


200o?/AT^ 


N(N -1)(N -4)(N -9) 


(6N^-27N^+42N^-30N+11) 


For n = (N+1 )/2 


2 2 
25o7Ar 


« tbo /Al » 2 

E[(L-y = ? V ? (3N^-18N +31) 


Acceleration error is given by 


E[(Jn-y„)'] = 


1440o?/aT^ 


N(N -l)(N^-4)(N^-9) 


[(2N+1)(9N+13) 


-70(N+l)n+70n‘^] 


For n = 1 and n = N 


E[(J„-y„)^] " 


1440o^/aT^ 


N(N‘^-1)(N -4)(N'^-9) 


(2N-1)(9N-13) 


12 



For n » (N+1 )/2 




EC(P„-y„)'] 


720Oq/AT^ 

N(N^-l)(N^-4) 


For .jerk we have 


ECCy-n-in) ] 


100 800o?/AT® 


N(N^-l)(N^-4)(N^-9) 


( 49 ) 


(50) 





13 


6. SUMMARY OF RESULTS FOR N LARGE 

Let n = N >> 1 . Then Table 1 summarizes the previously determined 
results. Note that T = (N-1 )aT=^ NaT. 

TABLE 1. SUMMARY OF RESULTS FOR n=N or n=l 


ORDER OF 
POLYNOMIAL 


I E[(y -y ) ] 

I — 




r 


E[(y„-y„) ] 


i^E[(jn-yn)^] 


0 

1 

4o^ 

~N^ 

2 

3 

2 

!g 

N 

9o^ 

__a 

N 

16o^ I 

N 


1920q/T^ 

1200Oq/T^ 

1 

1 

1 

N 

N 

ri 


2 4 

720oVT^ 

25920o^/T^ 

q 

! 

N 


j 

i 



100 BOOo^/T^ i 




N 1 




1 _ 


Let M = order of the polynomial fit. Then it appears that we may extrapolate 
our results to obtain the interesting equation for n = 1 and n = N. 


E[(y„-y„)^] = (H+D^o^/N 


(51) 


t 


1 


14 





Also we note from Table 1, that there is a large degradation in velocity accuracy 
when we go to the next higher order polynomial. Thus, when we desire a velocity 
estimate at the end point of the fit interval, we should use the minimum order 
polynomial fit. 

Now let n = (N+l)/2 » 1. Table 2 summarizes the previously determined 
results. 



(J 
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7. AN EXAMPLE 

Let us make up a hypothetical example in which we let 

y„ = w .in^ 


Using a random number generator, we can produce a table of noisy measurements, 
y* = yn + qn» as shown in Table 3. 


TABLE 3. LEAST-SQUARES FIT OP SECOND ORDER POLYNOMIAL 


n 

y 

q 

y * 

y 


y*-y 

y-y 

1 

1.1 

1.403 

2.503 

1.465 


1.038 

.365 

2 

1.4 

-.354 

1.046 

1.790 


-.744 

,390 

3 

1.9 

-.634 

1.266 

2.282 


- 1.016 

.382 


2.6 

.697 

3.297 

2.941 


.356 

.341 

5 

3.5 

-.686 

2.814 

3.763 


-.954 

.268 

6 

4.6 

1.375 

5.975 

4.762 


1.213 

.162 

7 

5.9 

.785 

6.685 

5.924 


.761 

.024 

8 

7.4 

-.963 

6.437 

7.254 


-.817 

-.146 

9 

9.1 

.759 

9.859 

8.751 


1 108 

-.349 

10 

11 

- 1.782 

9.218 

10.415 


- 1.197 

-.585 

11 

13.1 

-.600 

12.500 

12.247 

,_J 

.253 

-.853 


We will fit 


^n = 


a^ + a,n 
0 1 


+ a2n 


to the noisy y* data, using the equations of Section 4. 


) 


I 
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* 1 . 

4 


From Eq. 27 



I 


I 


I 


( 




1 

1 

1 

1 

1 

1 

1 

1 

1 

1 

1 


1 

2 

3 

4 

5 

6 

7 

8 

9 

10 

11 ' 


1 

4 

9 

16 

25 

36 

49 

64 

81 

100 

121 

Using the values 

of y* in 

Table 

3, we 

get 








61.600 ] 


P'^* =1 488.203 
4328.695 


Eq. 28 gives us the value of (P^P)"^ 

for N = 

11 of 



5174 

-1794 

130 ] 

1 

(P^P) ' = 

-1794 

759 

-60 


130 

-60 

5 


The estimated state vector is given by 



r . 1 

a^ 


1.3083 + 1.1 


0 



X = 

®1 

= (p^p)'VV = 

.0732 + .42 


^2 


.08375 + .034 


where the accuracy estimates were obtained by taking the square roots of the 
diagonal elements of the state error covariance matrix. 
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c = O^(P^P)'^ 

"q “ ' 

The res'Iduals are now estimated from 

y*-y = y* - 1.3083 - .0732n - .08375n^ 
n n n 

and their values are shown in Table 3. Note that 


11 


1 

n= 1 


(y* • y ) - -001 

'•^n ■'n' 


very commendable. 

If we have a large number of measurements, y^ will approach y^, making 

yt - y„ an estimate of q. Thus 
•'n n ^ 



a very good estimate for just 11 points, only off by 9% from the true value of 

A A ^ 

1. While Oq appears to be quite accurate, the individual values of q^ = y*-y^ 
are not nearly so accurate. It is desirable to have an equation to express 
the accuracy of estimating in the above manner, since it is frequently 
used to obtain the noise statistics for measurement data. Perhaps such a 
relationship will be found with the desirable accuracy. 

A A 

The last column of Table 3 shows y-y, the error in the estimate of y. 

As predicted by our previously derived equations, y is more accurate near the 
middle of the fit interval than it is near the end points of the interval. That 
is, from Eqs. 30 and 31 


18 




E[{y,|-y,,)^]= i .76 


Vt[(yi-yi)^]’ = 

\liUy(,-y(,)^i = i .46 

8. CONCLUSinNS 

One final thought to close this report. A lot of data are needed to 
reduce the measurement error standard deviation by a significant amount, 
say by a factor of 10. In this example with 11 points of data, we saw 
that at the end points, was reduced by a factor of .76, not much of 
a reduction. To get a factor of 10 reduction would have required about 
900 data points, a lot of points. 




$ 
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